
THE INCOMPLETE INVERSE AND 
ITS APPLICATIONS TO THE LINEAR 
LEAST SQUARES PROBLEM 

Georg E. Morduch ~ - 

Prepared by 

OLD DOMINION SYSTEMS, INC. 

Gaithersburg, Md. 

/or Goddard Space Flight Center 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • AUGUST 1977 


tech library kafb, nm 



TECH LIBRARY KAFB, NM 

miiiuiii 

□ Qtl.1,71 


1. Report No. 

NASA CR-2892 


2. Government Accession No. 


4. Titfe and Subtitle 

0?HE INCOMPLETE IWERSE ANT) ITS APPLICATIONS 
TO THE LINEAR LEAST SQUARES PROBLIM 


7. Author(s) 

Georg E. Morduch 

9. Performing Organization Name and Address 

Old Dominion Systems, Inc. 
Gaithersburg, Maryland 


12. Sponsoring Agency Name and Address 

National Aeronautics & Space Administration 
Washington, D. C. 20^46 

15. Supplementary Notes 


3. Recipient's Catalog No. 


5. Report Date 

August 1977 


6. Performing Organization Code 


8. Performing Organization Report No. 


10. Work Unit No. 


11. Contract or Grant No. 

NAS ^-21 '^87 


13. Type of Report and Period Covered 


14. Sponsoring Agency Code 


16. Abstract 

A modified matrix product is defined and it is shown that this 
product defines a group. The inverse of the group is called the in- 
complete inverse. Algorithms for computing the incomplete inverse are 
provided. 

The incomplete inverse has an important application to the least 
squares problem. It is shown that the incomplete inverse of an aug- 
mented normal matrix includes all the quantities (including the effect 
of "consider" parameters) associated with the least squares solution. 
In particular, an answer is provided to the problem that occurs when, 
(l) the data residuals are too large, (ll) there is insufficient data 
to justify augmenting the model by more than one term. A simple com- 
putation involving the incomplete inverse will tell which term will 
yield the best improvement in the data fit. This is of special inter- 
est to the processing of satellite data, where the model may always be 
augmented by any number of geopotential terms. The incomplete inverse 
[may thus be used to determine which geopotential terms most influence 
an orbit. 


17. Key Words (Selected by Author(s)) 

Least Squares 
Incomplete inverse 


18. Distribution Statement 


Unclassified - Unlimited 


Algorithms 

Matrix 


Cat. 67 

19. Security Classif. (of this report) 

20. Security Classif. 

(of this page) 

21 . No. of Pages 

22, Price 

Unclassified 

Unclassified 

8U 

$5.00 


*For sale by the National Technical Information Service, Springfield, Virginia 22161. 




ABSTRACT 


A modified matrix product is defined and it is shown that 
this product may be used to define a group. Each matrix in such a 
group possesses an inverse. To distinguish it from the regular ma- 
trix inverse, it is called the generalized inverse. The general so- 
lution of a matrix equation is expressed in terms of the incomplete 
inverse (a special form of the generalized inverse) . A minimum norm 
solution is derived, and it is shown that the associated matrix (the 
Penrose pseudo-inverse) is a generalized inverse. 

Algorithms for computing the incomplete inverse of a sym- 
metric matrix are derived. These algorithms, which are little more 
complicated than those for the regular inverse, utilize matrix sym- 
metry so that the matrix may be stored in upper triangular form. 
Similar use of matrix symmetry is made in the computation of the 
pseudo- inverse (the method thus requires only half as much computer 
core storage as the commonly used Andree algorithm) . 

The most important application of the incomplete inverse 
is to the least squares problem. It is shown that the incomplete 
inverse of an augmented normal matrix includes all the quantities 
(including the effect of 'consider' parameters) associated with the 
least squares solution. In particular, an answer is provided to the 
problem that occurs when, (i) the data residuals are too large, and 
(ii) the mathematical model may theoretically be augmented by a 
large number of terms, but (iii) there is insufficient data to jus- 
tify augmenting the model by more than one term. A simple computa- 
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tion involving the incomplete inverse will tell which term will 
yield the best improvement in the data fit. This is of special in- 
terest to the processing of satellite data, where the model may 
always be augmented by any number of geopotential terms. The in- 
complete inverse may thus be used to determine which geopotential 
terms most influence some observed orbit. 



PREFACE 


The incomplete inverse was first introduced in "An Exten- 
sion of NAP3.1F" (Morduch, 1975) in connection with the applica- 
tion of 'consider' parameters to the least squares solution. In 
that report it was defined in a rather complicated manner invol- 
ving matrix partitions and orthogonal transformations. For that 
reason the derivations and proofs of formulae involving the incom- 
plete inverse tended to be long and laborious , although the form- 
ulae themselves were relatively simple. The generalized inverse 
is defined in this report as the true inverse of a matrix in a 
group defined by a generalized matrix product, and includes the 
incomplete inverse as a special case. The reason for this report 
is not merely to give much simpler proofs for previously derived 
formulae, but also to describe the application of the incomplete 
inverse to the general and minimum norm solutions of matrix equa- 
tions . 


This effort was supported by NASA/GSFC under Contract No. NAS5-20532. 
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1 . 0 INTRODUCTION 


The reduction of observed data using the method of least 
squares is often fraught with difficulty. If this is caused by 
some unpredictable malfunction of a measurement instrument, then, 
provided that this is recognized, the offending data is excluded 
from the reduction. However, more frequently , the cause of the 
difficulty lies in the inadequacy of the mathematical modelling. 
This is often due not to any lack of knowledge of theory on the 
part of the investigator, but rather to the impracticability of 
having to solve for an extremely large set of unknown parameters. 
The concept of 'consider' parameters has been used by several in- 
vestigators to tackle the problem when this is so. 

The purpose of this report is to introduce a generalized 
inverse A’*’ of a matrix A, a particular form of which, the incom- 
plete inverse, has the property that its elements contain all the 
results of the least squares reduction including the effects of 
the 'consider' parameters. 

The generalized inverse is defined in Section 2 of this 
report, as is its particular form, the incomplete inverse. The 
general solution of a matrix equation is expressed in terms of the 
incomplete inverse. This is used to derive a minimum norm solu- 
tion. Although the method may inherently be less accurate than, 
say, the methods of Andree of Gram-Schmidt (Lefferts, 1969), all 
computations involve only symmetric matrices so that the computer 
storage requirement is much less than that of either of the two 


other methods. It is also shown that the matrix associated with 
the minimum norm solution belongs to the set of generalized inver- 
ses . 

An algorithm for the computation of the incomplete inverse 
of; a symmetric matrix is derived in Section 3. It is a generaliza- 
tion of the particular method of Gaussian elimination derived by 
the author for obtaining the usual inverse. (This method is used 
in the Navigation Analysis Program, NAPS. IF, used at Goddard Space 
Flight Center for orbit determination. Although the method involves 
exactly the same number of arithmetic computations as any other 
method of Gaussian elimination, which take advantage of the matrix 
symmetry, it is much simpler, and therefore easier to program, 
when the matrix to be inverted contains patterned zeroes.) 

Section 4 deals with the linear least squares problem. 

It is shown, amongst other things, that both the weighted least 
squares estimate, x, and the weighted sum of the squares of the 
residuals, C(x), are contained within the incomplete inverse, M', 
of a matrix M. The effect of consider parameters is taken into 
account in the derivation of a formula for the expected value of 

A A 

C(x). (C(x) is a measure of the goodness of the data fit and is 

therefore an important quantity to consider in the least squares 
reduction.) It is also shown that the change in the value of C(x) 
when a parameter is switched from a 'consider' to a 'solve' mode 
is given by a simple expression involving one multiplication and 
one ■ division. The last mentioned result is of considerable practi- 
cal significance for the following reason. Suppose that for some 
data reduction a given mathematical model results in an unaccept- 
ably bad data fit. Provided that the investigator knows how to 
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improve his model, he will then selectively (to augment the model 
by a large number of parameters in one step is not, in general 
desirable, since the results may be quite meaningless if too many 
parameters are solved for) augment his model by various parameters 
and then choose the one which gives the best fit. Such a task may 
be very arduous. However, by including as 'consider' parameters 
all such parameters as he might wish to include in his model, he 
may then select the one, which when switched from the 'consider' 
to the 'solve' mode yields the best improvement in the data fit. 

In other words , there is no need to perform a large number of data 
reductions using different mathematical models. One will do. 

The most important formulae are summarized in Section 5. 
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2.0 DEFINITION OF A 


GENERALIZED MATRIX INVERSE AND 
ITS APPLICATIONS TO THE SOLUTION 

OF MATRIX EOLATIONS 

- - * 


In this section we shall define a generalized matrix 
product between two square matrices. It will be shown that based 
on this generalized product a set of nxn matrices form a group. 

The generalized inverse of a matrix will then be defined as the in- 
verse of a matrix in this group. 

2. 1 The Generalized Matrix Product of Two Square Matrices. 

Definition : The generalized matrix product of 

any two square matrices (of equal dimen- 
sions) A and B is denoted by A-B and is 
defined by 

A-B = -ARB +A(I-R) +(I-R)B, (2.1) 

where I is the identity matrix and R is a square matrix satisfying 

R = R^ (2.2) 

and 

RR = R (2.3) 

It can easily be shown that for any three n^n matrices 

A, B and C 

(A-B) . C = A- (B-C) (2.4) 

k 



so that the generalized product is associative . 
It can also be shown that 


A- (-R) = (-R)-A = A, 


(2.5) 


for any n«n matrix A. Hence (-R) is the identity element . 


The inverse of A, which we shall denote by A , must satis- 


fy 


A-A* = -R 


( 2 . 6 ) 


As is well known from group theory, it follows from equa- 
tion (2.6) that 


A*-A = -R 


(2.7) 


It follows from the above that based on the generalized 
product, all n«n matrices that possess an inverse form a group. 

We hence conclude that the inverse is unique and also that 


(A*)* = A, 


and 


(A-B)* = B*-A* 


( 2 . 8 ) 


(2.9) 


It follows from equations (2.1) and (2.2) that 
(A-B)" = -B^RA^ + (I-R)A^ +B^(I-R), 


1 . e. , 


(A-B)^ = B^ -A^ 


( 2 . 10 ) 


Hence taking the transpose of equation (2.5) we deduce that 
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(A*)^ -A^ 


-R 


Since by definition 

(a')*. a' = -R, 

we conclude that 

(A*/ = (A^)* (2.11) 


It follows from the above equation that if A is symmetric 
then so is A*. 


2 , 2 The Generalized Inverse . 

A* has been defined as the inverse of A with respect 
to the generalized product based on R. Henceforth it will 
be referred to simply as the generalized inverse of A, However, 
when any risk of confusion arises, it will be referred to as the gen- 
eralized inverse of A with respect to R and will be denoted by A,R. 

We shall next derive some useful formulae involving A*. 

At this point it is convenient to define the square matrix S by 

S = I-R (2.12) 


It follows from the above definition and from equations 
(2.2) and (2.3) that 


and 


S = s', (2.13) 

SR = RS = 0, (2.14) 

SS = S (2.15) 
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Equation (2.1) may then be rewritten in the form 

A-B = -ARB + AS + SB (2.16) 

From the above and equation (2.6) we deduce that 

ARA* = R + AS + SA*, (2.17) 

whence 

A* = (AR -S)"^ (R +'AS) (2.18) 

Note that it follows from the uniqueness of A* that if A* 
exists then (AR-S) must be non-singular. From equations (2.16) and 
(2.7) we similarly find that 

A*RA = R + A*S + SA, (2.19) 

whence 

A* = (R + SA)(RA -S)"^ (2.20) 

It follows from equations (2.8) and (2.18) that 

A = (A*R -S)"^ (R + A*S) , (2.21) 

and similarly from equations (2.8) and (2.20) that 

A = (R + SA*) (RA*-S)'^ (2.22) 

Pos tmul.tiplying equation (2.17) by R we obtain 

ARA*R = R + SA*R, (2.23) 

whence 

(AR-S)A*R = R (2.24) 


= I, we therefore conclude 

T 


Since (AR-S)S 


-S and R+S 


that 


(AR-S) (A*R-S) = I 


Hence also, 


(A*R-S) (AR-S) = I 


Pr emu It ip lying equation (2.23) by R we find that 

(RAR)(RA*R) = R, 


whence 

(RAR+S) (RA*R+S) = I 


Since 


AR-S = A-(A+I)S, 

we deduce from equation (2.25) that 

A(A*R-S) = I + (A+I) (SA*R-S) 


(2.25) 


(2.26) 


(2.27) 


(2.28) 


(2.29) 


•k 

2.3 The Rank of A when SA S = 0. 


In general, if M, B and C are three nxn matrices satisfy- 


ing 


and 


M = B + C, 


BC = 0, 


(2.30) 


(2.31) 


then clear).y 
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rank(M) < rank(B) + rank(C) 


(2.32) 



nncl 


rank(B) < n-rank(C) 

If furthermore, M is non-singular so that rank(M) ■ 
then it follows from (2.32) that 

rank(B) > n-rank(C) 

Comparing inequalities (2.33) and (2.34) we deduce 
rank(B) = n-rank(C) 


Since 


I = R+S and kS = 0, 

we therefore conclude that 

rank(S) = n-rank(R) 

Since for any matrix M, 

(I-RMS) (I+RMS) = I, 

it follows that I-RMS is non-singular. 

Since 


I-RMS = S-RMS + R 


we deduce that for any matrix M 

rank (S-RMS) = n-rank(R) 


(2.33) 
n, 

(2.34) 
that 

(2.35) 


(2.36) 

(2.37) 


(2.38) 


It follows from equation (2.28) that 
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RAR + S is non- singular . 


Since 

RAR + S = (S-RAS) + RA, 

we hence deduce that 

rank(RA) = n- rank (S-RAS) 

Since equation (2.38) is valid for any matrix M it follows 
from the above that 

rank(RA) = rank(R) (2.39) 

Since quite obviously 

rank(RA) < rank (A) , 

we conclude that 

rank(R) < rank (A) (2.40) 

Similarly 

rank(R) < rank (A*) (2.41) 

We shall now make use of the assumed relation 

SA*S = 0 (2.42) 

Pos tmul tip lying equation (2.17) by S, we deduce with the 


aid of the above equation that 



Hence 


rank (A) < n-rank(S-RA*S) 

Letting M = in equation (2.38) we obtain 
rank(S-RA^^S) = n-rank(R) 


(2.44) 


(2.45) 


Hence 


rank (A) < rank(R) (2.46) 

Comparing inequalities (2.40) and (2.46) we conclude that 

rank (A) = rank(R) (2,47) 

Note that inequalities (2.40), (2.41) and (2.45) follow 
from the existence of A*, whereas inequality (2.46) is a conse- 
quence of equation (2.42). 


2 . 4 An Orthogonal Transformation of the Generalized Inverse . 

Let V be an orthogonal matrix satisfying 


w' = vW = I (2.48) 

Let Aq , Rq and So be defined by 

A^ = vAv^ , (2.49) 

Ro = vRv^ , (2 . 50) 

and So = vSvl (2.51) 
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It follows from the above, that 


and 


K = Ro‘ 
R„R„ = Ro 
S.. = I-R„ 


Writing equation (2.17) in the form 


(2.52) 

(2.53) 

(2.54) 


(2.55) 

We deduce after premultiplying the equation by v and 
postmutiplying it by , and inserting v'^v(=I) between all ma- 
trices that 


(A„R„-S, )v(A*R)v^= R„+A„S, , 

whence 

v(AtR)v'= (A„Ro-S„ )"' (R„ + AoSo) (2.56) 

Since the right hand side of the above equation equals 
Ao ,Ro by definition, it follows that 


aJ.Ro = v(AtR)v^ 


(2.57) 


2 . 5 The Incomplete Inverse . 

The incomplete inverse is defined as the generalized in- 
verse with respect to R, where R is a diagonal matrix. 

It thus follows from equation (2.3) that the elements of 
R are either zero or one. The elements of S are similarly either 
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zero or one. 


The reason for calling the inverse incomplete is as fol- 
lows. Inverting equation (2.21) we obtain 

A“' (R + A*S)“' (A*R-S) (2.58) 


The inverse of A may thus always be obtained directly from 
the generalized inverse. However, the total nuinber of computations 
in obtaining A“ ^ is independent of whether or not incomplete in- 
verses are obtained in intermediate steps. 


Let us now assume that we are given matrices A and R, 
where R is a diagonal matrix, whose elements are either zero or 
one . 


We can then find an orthogonal matrix v satisfying equa- 
tion (2.48) such that R„ . as given by equation (2.50), may be writ 
ten in piirtitioncd form as 


Ro 


0 0 

0 I 


(2.59) 


Correspondingly let A„ and So (as defined by equations 
(2.49) and (2.51)) be written as 


and 



(2.60) 



0 

0 


(2.61) 
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It follows from the above that 


and 


(A,R,-S„) = 


-I 

0 


H 

L 


(2.62) 


(R, + = 


G 

K 


0 

I 


(2.63) 


Hence , 


Since 


(AoRo-So ) 


- 1 


•I 

0 


HL 
_ - 1 


(2.64) 


A*,R, 


= (A,Ro-So)’'(Ro + A„S,), 


we obtain, when substituting the expressions from equations (2.63) 
and (2.64) on the right hand side of the above equation 




a: . R„ = 


(HL^K-G) HL" ^ 


L"^K 


- 1 


Since by equations (2.48) and (2.57) 


(2.65) 


A*,R = vUAr,Ro)v. 


( 2 . 66 ) 


it follows from equation (2.65) that 


A^R = 


(H"^LK-G) HL"^ 


L"^K 


- 1 


V 


(2.67) 



where 


G H 

K L 


vAv^ 


( 2 . 68 ) 


It can be seen from the above that if A is appropriately 
defined then A* will yield matrices of the form L”' , L~' x and 
L~' B even if in a jumbled (to the extent that v differs from the 
identity matrix) form. The application of the incomplete inverse 
to the solution of matrix equations involving consider parameters 
will be dealt with in a later section. 


2 . 6 The General Solution of the Matrix Equation Ax = y , 

Where A is a Symmetric Matrix . 

Let us consider the matrix equation 

Ax = y, (2.69) 

where A is a symmetric nxn matrix and x and y are n-vectors. 

Let us assiame that rank (A) = r. Let R be defined as a diagonal ma- 
trix, whose elements are either zero or one. We can clearly choose 
R such that the non- zero columns of AR are linearly independent 
and the remaining columns of A are linear functions of the columns 
of AR. In other words, there exists an nxn matrix P such that 

AS = ARP, (2.70) 

and 

SP = 0 , (2.71) 
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where 


S = I-R. 


Since the non- zero columns of AR are linearly independent 
it follows that if there exists a vector u such that 



ARu 

= 0 , then Ru = 0 . 

(2.72) 


We shall now show that AR-S is non- singular . 

Let us 

assume 

that there exists 

a vector u such that 




(AR-S)u = 0 

(2.73) 


Premultiplying 

the above equation by R and 

S, respec- 

tively 

we deduce that 





RARu = 0 

(2.74) 

and 






(SAR-S)u = 0 

(2.75) 


Pos tmultiplying the transpose of equation (2.70) by Ru 

we obtain 


SARu = P^RARu, 
whence by equation (2.74), 

SARu = 0 (2.76) 

Adding equations (2.74) and (2.76) we obtain 

ARu = 0. 
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It hence follows from equation (2,72) that 


Ru = 0 (2.77) 

Adding equations (2.76) and (2.77) and subtracting equa- 
tion (2.75) from the result we find that 

u = 0 (2.78) 

We have just shown that (AR-S)u = 0, implies that u = 0. 
Consequently (AR-S) is non-singular. It hence follows from equation 
(2.18) that A* exists. 

Subtracting equation (2.71) from (2.70) vie obtain 

(AR-S)P = AS, (2.79) 

whence 

P = (AR-S)"' AS (2.80) 

The right hand side of the above equation equals AS, 
as may be readily verified by postmultiplying equation (2.18) by 
S . Hence 

P = A*S ■ (2.81) 

It then follows from equation (2.71) that 

SA*S = 0 (2.82) 

Hence equation (2.43) applies, i.e., 

A(S-RA*S) = 0 (2.43) 

Taking the transpose of the above equation we obtain 

IT 



0 


(2.83) 


(S-SA*R)A = 

Therefore, in order that equation (2,69) have any solu- 
tion at all, it is necessary that 

(S-SA*R)y = 0 (2.84) 

From the above and equation (2.29) it then follov7S that 

A(A*R-S)y = y (2.85) 

Comparing equations (2.69) and (2.85) we conclude that 

= (A*R-S)y (2.86) 

is a solution of equation (2.69). In view- of equation (2.43) the 
general solution of equation (2.86) may be written as 

X = Xo + (S-RA*S)a, (2.87) 

where the vector a satisfies 

Ra = 0, (2.88) 

, but is otherwise arbitrary. 

Adding equations (2.84) and (2.86) we obtain the following 
alternative form for x^ . 

X, = RA*Ry (2.86)'. 

That X as given by equation (2.87) is the most general 

solution of equation (2.69) can be seen from the follovzing argum.ent. 

Since by assumption rank (A) = r, the general solution of equation 
l8 ' 



(i*. .69) muMl. conLain an arbitrary linear combination of (n-r) linear- 
ly independent vectors. Since by equation (2.47) rank(R) = rank(A) 
it follows that rank(R) = r, whence by equation (2.45) 
rank(S-RA*S) = n-r. x therefore is of the required form and is 
hence the general solution. 

Note that we have shown that 

(1) The general solution of equation (2.69) is 
given by equation (2.87), where x^ is any 
solution of equation (2.69), 

(ii) Given that SA*S = 0 it follows that 

(A*R^S)y is a solution of equation (2.69) 

and (iii) There exists an R and a corresonding S, 

both of which are diagonal matrices, and 
which satisfy SA*S = 0. 


It will later be shown (see Section 3, equations (3.10)', 
(3.12)', and (3.13)') that if we define the symmetric matrices M and 
Rm by 



(2.89) 


(2.90) 



(2.91) 
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where x^, is given by equation (2.86) or equivalently be equation 

(2.86) '. Since it follows from equation (2.86)' that 

Rx„ = x„ (2.92) 

equation (2.91) may also be written in the form 

(2.93) 

The significance of the term (y^x^-C) will be explained 
later in Section 4 dealing with least squares parameter estimates. 

Finally we note that it follows from equations (2.84) and 

(2.87) that 


M 


y^x„-c 


y^x = y^x„ 


(2.94) 


In other words, the product y’^ x is the same for all 
solutions X of equation (2.69) 


2.6.1 The Minimum Norm Solution 

If A in equation (2.69) is singular then S is not equal 
to zero and equation (2.87) will yield an infinite number of solu- 
tions for equation (2.69). Of all possible solutions, the minimum 
norm solution is the one that minimizes the norm Nq(x) given by 

Nq(x) = x' Q” X, (2.95) 

where Q is a diagonal positive-definite matrix, which defines the 
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norm. Let us denote the value of x that minimizes Nq(x) by Xq . 
Since x is a function of a it follows that the parital deriva- 
I ivc's ol' M(j(x) with respect to a must vanish. We hence find with 
I he aid o( (‘(((i.il ion (2.87) I hat 


(S-RA*S)' Q'^ [ x„ + (S-RA*S)a] = 0, 

i.e. , 

(S-SA*R)Q"^ [ x„ + (S-RA*S)a] = 0 (2.96) 


Since S, R and Q~^ are diagonal matrices, they commute 
and the above equation may be simplified to 


(SA*RQ"*RA*S + Q“^S)a = - (S-SA*R)Q‘^ x^ (2.97) 

Since it follows from equation (2.88) that 

Sa = a (2.98) 

equation (2.97) may be rewritten as 

(SA*RQ"^ RA*S + Q"^ )a = - (S-SA*R)Q“^ x„ (2.99) 


fine 


In order to simplify the algebra it is convenient to de- 

H = RA*S (2.100) 


It follows from the definition and the properties of R 
and S that 


and 


HH = SH = HR = 0 , 

RH = HS = H 


( 2 . 101 ) 

( 2 . 102 ) 
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Equation (2.99) may now be re-expressed in the form 


(H"Q~'H + Q')a = -(S-hMQ“'x„. (2.103) 

(H^Q + Q being the sijm of a positive-definite and a semi- 
positive-definite matrix must be non-singular. Hence 

a = -(H^Q'^H + Q"')“' (S-HMQ"'x„ (2.104) 

From the above and equations (2.87) and (2.100) it follows 

that 

Xq = NXo , (2.105) 

where 

N = I -(S-H) (h"q~^H + Q~^)"^ (S-H^)O'^ (2.106) 

An alternative form for N is obtained as follows. By 
definition of the inverse 

(H^Q“^ H + )'^ (H^Q"^ H + Q“^ ) = I (2.107) 

Pos tmultiplying the above equation by Qtf we obtain 


(h"q“^ H + Q“^ ) ^ H^Q"^ (HQH" + Q) = Qh\ (2.108) 

whence 

B = QH^ (HQH^ + Q)‘^ (2.109) 

where 

B = (H^ Q"^ H + Q"^ )“^ Q"^ (2.110) 
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From equations (2.107) and (2.110) we deduce that 

(H^ Q‘^H + Q"^ )‘^ = Q -BHQ (2.111) 



It follows from equation (2.109) that 


HB = I -Q(HQH' + Q)’^ (2.112) 

Writing equation (2.106) in the form 

N = I -S(irq"^H + Q“^ )“^ SQ"^ + (S-H)(H^Q"^H + Q“^ )H^ Q”^ 
+H(H^Q’^H + Q‘^)“^ SQ"\ 

we deduce with the aid of equations (2,110) and (2.111) that 

N = I -S(Q-BHQ)SQ“^ + (S-H)B + py SQ"^ (2.113) 


Since S and commute it follows from the above and 

equations (2.109) and (2.112) that 


N = I-S + BH + B-HB + QB^ Q“^ 

= -S + QH" (HQH^ + Q)“^ (I+H) + Q(HQH' + Q)"^ 

+ Q(HQH^ 4- Q)‘^ H 

Hence 


N = -S + (Q + QH^)(HQH^ + Q)"^ (I + H) 


(2.114) 


Note that since it follows from equation (2.86)' that 
Sxo = 0, whence also HXo = 0, equations (2.105), (2.106), and 

(2.114) may be replaced by 

\ 

Xq = x„ + (S-H)(h"q"^H + Q“^)'^ H^Q"^x„ (2.115) 

and the alternative form 

Xq = (Q + oy )(H0H" + Q)'^ Xo (2.116) 
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The choice as to which equation to use may be made by 
noting that + Q”^) is effectively an (n-r) x (n-r) matrix and 

(HQH^ + Q) is effectively an rxr matrix. (Both matrices are n*n 
matrices and r is the rank of A) . 


2.6.2 Relationship of the Minimtim Norm Solution to the General- 

ized Inverse . 

It follows from equation (2.110) that BS = 0. We there- 
fore deduce from equation (2.112) that 

(HOHW Q)"^ S = Q"^ S (2.117) 

Since (I+H) (S-H) = S it follows from the above and equa- 
tion (2.114) that 


N(S-H) = -S + (Q+QH')Q'^ S, 


i.e. , 


N(S-H) = 0 


(2.118) 


From the above and equation (2.106) we deduce that 


NN = N 


(2.119) 


Since be equations (2.43) and (2.100) 


A(S-H) = 0, 


( 2 . 120 ) 


we similarly deduce from equation (2.106) that 


2k 


AN 


A 


( 2 . 121 ) 



Since by equations (2.84) and (2.100) 


(S-H")y = 0, (2.122) 

we also deduce that 

N"y = y (2.123) 

Let us now define Qo as the diagonal matrix, whose ele- 
ments equal the square-root of the corresponding elements of Q, so 
that 


QoQo = Q 


(2.124) 


Also define 


Ro = q;' NQ, 


(2.125) 


It can be seen from equation (2.106) that Ro is symmetric, 


1 . e. , 


Ro = Rc 


(2.126) 


It can easily be seen from the definition of Ro and from 
equation (2.119) that 


Rq " Rc 


(2.127) 


Defining So by 


So I -R o , 


(2.128) 


it follows that 


• So So = Sc 


(2.129) 
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and 


SoRo = RoSo = 0 (2.130) 

Let L be defined by 

L = R,q;'rA*RQ:^Ro, (2.131) 

and let us consider the expression 

E = (QoAQoRo -So)(LRo -So) 

It follows from equations (2.125) through (2.132) that 

E = QoAQoL + So . (2.132) 

= QoANRA*RQr^ Ro + So 

From the above and equations (2.121) and (2.23) we de- 
duce that 

E = Qo (R + SA*R)Qo"^Ro + S„ 

= Ro + So + Qo (H^ -S)Qo"' Ro 
= I + [RoQ:' (H-S)QoJ 
= I + [q:*N(H-S)QoJ' 

= I , by equation (2.118) 

It follows from the above and equation (2.132) that 

(QoAQoRo -So)"^ = LRo -So (2.133) 

From equations (2.121) and (2.125) we deduce that 
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(QoAQo)R 


Qo AQo ) 


(2.134) 


whence 

(QoAQ.,)So = 0 (2.135) 

From the above and equation (2.18) we conclude that the 
generalized inverse of (Qo AQo ) with respect to Ro is given by 

(QoAQo)tRo = (LRo -So)Ro, 

whence by equations (2.127), (2.130), and (2.131) 

(QoAQo)tRo = L (2.136) 

We shall next show that Xq = (QoLQo)y. 

It follows from equations (2.123) and (2.125) that 

Qo"'RoQoy = y (2.137) 

From the above and equation (2.131) we deduce that 

LQoy = RoQ^'RA^Ry 

= RoQ:'x„, (2.138) 

by equation (2.86)' 

From the above and equation (2.125) we find that 

V LQoy = Q^^NXo 

= Qo Xq 

by equation (2.105) 

Hence 

Xq = (QoLQo)y, (2.139) 
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where L is given by equation (2,136) as the generalized inverse of 
QoAQo with respect to Ro. In other words we have shown that there 
exists a generalized inverse, which is directly related to the mini- 
mum norm solution Xq . 

It is interesting to note that L is the Penrose pseudo- 
inverse of (QoAQo). (See e.g., (Lefferts, 1969)). 

Since it was shown that E in equation (2.132) equals the 
identity matrix it follows that 

(QoAQo )L = Ro (2.140) 

Since Ro is S3?mmetric we deduce that 

(QoAQo )L = L(QoAQo) (2.141) 

From equation (2.131) we find that 

LRo = L, 

whence we deduce from equation (2.140) that 

L(QoAQo)L = L (2.142) 

Similarly we deduce from the transpose of equation 
(2.134) and from equation (2.140) that 

(QoAOo)L (QoAQo) = QoAQo (2.143) 

Since equations (2.141) through (2.143) are the equations 
satisfied by the Penrose pseudo -inverse of the symmetric matrix 
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(QoAQo) it follows that L is the Penrose pseudo -inverse of (QoAQ^) , 


2 . 7 A Relationship Between Two Incomplete Inverses . 

Let U and V be defined by 

U = A*R , (2.14A) 

and 

V = A^R' , (2.145) 

where R and R' are diagonal matrices, whose elements are either 
zero or one. Let us assume that R' has more noh-zero elements than 
R, so that we may write 


R’ = R + K. 

where the elements of K are either zero or one, and 


(2.146) 


RK = 0 


(2.147) 


It follows from equations (2.21) and (2.18), respectively. 


that 


A = (UR -S)“^ (R + US) , 


and 


V = (AR* -S')‘^ (R + AS') 


(2.148) 


(2.149) 
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From equation (2.148) we deduce that 


and 


AR’-S' 


AS' + R’ 


= (UR-S)"* |^(R+US)R' -(UR-S)S'J 

= (UR-S)"^ [(R + US)S' + (UR-S)R'] 


(2.150) 


(2.151) 


It follows from the above and equation (2.149) that 
V = [ U(SR'-RS’) + (RR' + SS') U(RR' + SS ' ) - (SR ' -RS ' ) ] (2.152) 


Since 


S' = I-R' 


we deduce from equation (2.146) that 


S' = S-K 


(2.153) 


It follows from equation (2.147) that 


SK = K 


(2.154) 


We hence deduce that 

RR' + SS' = R + S-K 


= I-K, 


and SR'-RS' = K, 

whence 


V = UK + (I-K) [^U(I-K)-kJ . (2.155) 
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The above equation gives the desired relationship between 
the two incomplete inverses defined by equations (2.144) and (2.145) 

Note that if K = S, so that R' = I and hence V = A"^ , then 

V = (US + R)"^ (UR-S) (2.155) 

The above equation could also have been obtained directly 
by inverting both sides of equation (2.148). 


31 



3.0 AN ALGORITHM FOR COMPUTING 
THE INCOMPLETE INVERSE 


Let M, M*, R, and S be partitioned according to 


M 


C 

Y A 


M 


D 

X B 


(3.1) 


R = 




S 



where M and M are symmetric matrices and R and 
matrices, whose elements equal either one or zero, 


0 


(3.2) 


S.J , 

S are diagonal 
and R + S = I. 


Since by definition the incomplete inverse M must satis- 
fy equation (2.17) it follows that 


1 

n 

—1 1 


1 

n 

o 

1 


D X'' 

1 

> 

1 


0 R„ 


X B 



1 

A 

o 
1 


c y " 


Se 0 


Se 0 


D X^ 

= 


+ 









1 

o 

a 


1 

< 

>• 

1 


1 

o 

CO 

o 
1 


0 Se. 


X B 


Equating the partitions we hence obtain 

CR, D + Y^R„X = Re + CSc + Se D, (3.4) 

YReD + ARoX = YSc +SaX, (3.5) 
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YR,X" + AR„B 


R„ + AS, + S, B 


(3.6) 


Note that since matrices M, M‘, R, and S are symmetric, 
the equation for the fourth partition of equation (3.3) is redundant. 

Since by equation (2.25) 

(AR,-S,)-' = A*R,-So, (3.7) 


we deduce from equation (3.5) that 

X = (aX-So)Y(S, -R,D) 


(3.8) 


In order to simplify the algebra let us write 


F = -(A*R,-S, )Y, (3.9) 

so that 

X = F(R,D-Sc) (3.10) 

Substituting the expression for X from equation (3.10) on 
the left hand side of equation (3.4) we obtain 


(C + y'r,F) (R,D-Sc ) = Re + Se D (3.11) 

Comparing equations (3.11) and (2.20) we see that 

C + Y^R,F = D* 

We hence deduce with the aid of equation (2.8) that 

D = (C + y'r,F)* (3.12) 
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It follows from equation (3.6) that 


(AR„-S,)B = R, + AS, -YR, X', 

whence by equations (2.18) and (3.7) 


B = A* -(A*Ra-S,)YR,x' 



From the above and 

equation (3.9) we finally obtain 



B = 

^ A* + FR,X^ 

(3.13) 


Note that if C, D, Rc 

and are scalars, X and Y vectors and 


Rc = 

0 and Sc ■ = 1 


then 

equations (3.9), (3.10), 

(3.12), and (3.13) reduce to 



X = 

(A*R, -S, )Y , 

(3.10) ' 


D 

= Y^R,X -C 

(3.12) ’ 

and 

B 

II 

> 

(3.13) ' 

3.1 

The Computation of 

the Incomplete Inverse 



Let us define the symmetric matrix M by 


M 


M' 

F 



(3.14) 


= C + Y^R.F 


where 
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M' 


(3.15) 



It follows from the above and equations (3.10), (3.12) 
and (3.13) that the partitions of 


M* 


D 

X B , 


are given by, 



D = (M')* 

(3.16) 


x' = (DRc -S, )f' 

(3.17) 

and 

B = A* + FR, X^ 

(3.13) 


It can thus be seen that M* may be computed directly from 
M. (M')* may be computed by partitioning M' similarly to the way M 
was partitioned. This process may be incorporated in a procedure as 
follows . 


If M is a sjmimetric nxn matrix, define M„ by 


Then for r = n, n-1, 
rxr matrix according to 


= M (3.18) 

n-2,...2 partition each S3nnmetric 


and 


M, 


Also define 



A, 


R = R, 


S" = S 



K 0 


s; 0 

II 

r 

, = 

r 


o 

0 


0 So 


(3.19) 


(3.20) 


(3.21) 
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Corresponding to equations (3.9) and (3.15) define 

f; = -y7(R:a*-S;) (3.22) 

M,_, = C, +y;r:f„ (3.23) 

Also define 

R^"' = r; and S'~' = S/ (3.24) 

and A, by 

A, = M, , (3.25) 

whence 

= A* (3.26) 

Equations (3.18) through (3.26) are referred to as the 
elimination equations. The inversion of M is completed through the 
use of back- substitution equations (3.27) through (3.29). 


For 


r = 2, 3, 


n 


Mr 


M’ 


r - 1 

Xr 


x; 

Br 


where corresponding to equations (3.17) and (3.13) 


Xr' = (M*j r; -s;)f/ , 


and 


Br 


(3.27) 


(3.28) 
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A* + Fr R; Xr' 


(3.29) 



Finally, since M„ = M it follows that 

M* = M* (3.30) 

NOTE 1 

If M is calculated on a digital computer, then since all 
matrices M, are symmetric, only the upper triangular portion of each 
matrix need be stored in computer memory. Also, storage space may be 
shared by the following matrices : 

(y/, f/. xJ), (A,, A*, B, ) and (C, , 

In view of equation (3.28) n temporary storage locations must how- 
ever be allocated. The total number of storage locations required to 
compute M* from M is thus n(n+3)/2. 

NOTE 2 

Since is a 1»1 matrix and the elements of R are 

either zero or one, equations (3.22) and (3.23) may be simplified for 
the two cases Rj = 1 and Rj = 0. Similarly A* as given by equa- 
tion (2.18) may be simplified. Corresponding to the two cases we 
obtain: 


If RJ = 1 then 

(3.31W) 
(3.22W) 
(3.23W) 


and if R^' = 0 then 

A* = -A, , 



(3.31V) 
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(3.22V) 


and 




= c, 


(3.23V) 


NOTE 3 

Although A, and Br in equations (3.19) and (3.27) are 
assumed to be 1»1 matrices, equations (3.18) through (3.30) remain 
valid if e.g., they are kxk matrices. This fact may be utilized, 
when the matrix M is so large, that only part of it may be stored in 
main computer memory. 

NOTE 4 

If, e.g.. 



\ 

o 

M 

1 


— 1 
o 

1 

M = 


and R = 



1 0 


1 

o 


then M* = M~' = M. However, the above inversion scheme will ob- 

viously break down. It will next be shown that the inversion scheme 
can always be applied when the matrix (RMR+S) is positive-definite. 

3 . 2 A Sufficient Condition for the Applicability of the 

Inversion Scheme . 

It will now be shown that if the matrix (RMR+S) is posi- 
tive-definite, then the inversion scheme is always applicable. 

First we shall show that M* exists if and only if (RMR+S) is non- 
singular . 

It may readily be verified that 
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(R-S + SMR)(R-S + Sm) = I, 

SO that (R-S + SMR) is non-singular. Since 

(MR-S) = (RMR + S)(R-S + SMR), 

it follows that (MR-S) is non-singular if (RMR+S) is. We hence, 
deduce with the aid of equation (2.18) that M* exists if (RMR+S) 
is non-singular. Since it follows from equation (2.28) that 
(RMR+S) is non-singular if M* exists, we conclude that M* exists, 
if and only if (RMR+S) is non-singular. 

To show that the inversion scheme always can be applied 
when (RMR+S) is positive-definite, it is, in view of the recursive 
nature of the scheme, sufficient to show that, when M is parti- 
tioned according to equation (3.1), both (R„ARa+Sa) and (RcM'R^+S,. ) 
l^where M' is given by equation (3.15) j are positive-definite. 

It follows from equations (3.1) and (3.2) that 

R, CRc +Sc Rc Ra 

RMR+S = (3.32) 

RaYR, R,ARa+S„ 

Since RMR+S is positive-definite, it follows that 


Ro ARo + Sa 


= R, CRc + Sc - RcY^R^ (R„ARc + S„ ) "' R^YRc (3.33) 


are positive-definite. jThat G must be positive-definite can be 
seen from the following. Let the vector v be defined by 
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-u^ R (RoARa + Sa , whence (RMR + S) v = u^Gu. 

The result thus followsj . 

An alternative expression for G is obtained with the aid 
of equation (2.28) as 

G = R, CR, + S, -R, Y" R„ (A*R, 4-S, )R„YR, , 

i . e . , 

G = R, [C -Y^R„A*R,Y] Rc + Sc 

= Rc [C -Y'R„(A^'^Ra - S.,)Y]Rc + Sc 
= Rc [C + Y' R„F] Rc + Sc , 

by equation (3.9) 

It follows from the above and equation (3.15) that 
G = RcM'Rc + Sc . We have hence shown that both RoARa + So and 
RcM'Rc + Sc are positive-definite and that therefore, the inversion 
procedure always can be applied when RMR+S is positive-definite. 

3 . 3 Computation of the Determinant of the Normalized Inverse . 

In this sub- section we shall compute the determinant of 

M = A~^ (RMR + S)~^ a"\ (3.34) 

where A is the diagonal matrix, whose non-zero elements equal the 
square-root of the corresponding elements of (RMR+S)" ^ . It follows 
from equation (3.34) that 

det(M) = det(RMR + S)"^ /(det a)^ , 

ho 



and hence that 


det(M) = l/det(RMR + S)det( A^) (3.35) 

It follows from the definition of A, that A^ is the di- 
agonal matrix, whose non-zero elements equal the corresponding diagon- 
al elements of (RMR+S)~^. Since by equation (2.28), (RMR+S)"^ = 

(RM*R+S) , it follows that 


det(M) = det(RM*R + S)/det( A^ , 


(3.36) 


where det( A^ equals the product of the diagonal elements 


of (RM'"R + S) 


It follows from equations (3.1) and (3.2) that 


RM*R + S = 


whence by equation, (3.13) 


RM*R + S = 


R, DR, + Sc 
R.X Rc 


+ 


RcDRc + Sc 
RoX Rc 


0 

0 


RcX' R„ 

Ra BRo + So 


0 

R„A*Ra + S„ 

Rc X^ Ro 
Rc FRc X^ Rc 


(3.37) 


(.3.38) 


Since X = F(RcD-S, ) it follows that 


Re X' Rc = Rc DRc F' Rc , 


1+1 



and 

whence 


R„FR, X R„ = 


RcXRc 


R„FR, DR, F' R„ , 
Rq FR, DRc , 


R, X' Ra 


RcDR, + S, 

R, FR, X' R„ 


R.XR, 


Rr F R. 


(3.39) 


Since the value of a determinant is unaffected by the addi- 
tion to any one column of a linear combination of the remaining 
columns, we deduce from equations (3.38) and (3.39) that 




det(RM R + S) = det 


RcDR, + S, 
Ro XRc 


0 


R„A*R, + S„ 


Hence 




det(RM R + S) = det (R, DR, + Sc)det(R«A R„ + S, ) . 
From the above and equation (3.16) we deduce that 


det(RM*R + S) = det (R, (M' )*R, + S, )det(R„A*Ra + S, ) (3.40) 


Applying equation (3.40) to matrix M, as defined by equa- 
tion (3.19) we obtain with the aid of equations (3.21) through (3.24) 


det(R'M*R' + S' ) = det(R"^ R'"^ -f S'"^ ) 

X det(Rl A* R',, + Sa ) 


(3.41) 
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Since M„ = M we deduce from the above and equation (3.26) 


that 


det(RM*R + S) = IT det(R;A*R;. + S^), (3.42) 


where 


r\ = R^ and Sa = 


Since Ro , and are 
simplified to 


1x1 matrices, the above may be 


det(RM*R + S) = ~1 T(RoA*R'<. 4- SL ) (3.43) 

r=l 

2 

It follows from equations (3.26) and (3.27) that det(A ), 
which equals the product o£ the diagonal elements of (RM R + S) , is 
given by 


det(A") = (rU*Ro + S,') it (R^rR'a + SL) (3.44) 

r*2 

It hence follows from equations (3.36), (3.43), and (3.44) 
and, since R^ equals either one or zero, that 


where 


det(M) = TT Qr , 


A*/B, if Ri = 1 
1 if rI = 0 


(3.45) 


(3.46) 


Note that, since for large matrices det(M) may be extremely 
small, det(M) is best computed logarithmically. 
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3.4 


The Determination of R when M is a Semi-positive-definite 


Matrix . 

In this subsection it will be shown how a diagonal matrix 
R is determined such that the non- zero columns of MR are linearly 
independent, but the remaining columns of M are linear functions of 
the columns of MR. 

It was previously shown (see Section 2.6) that if R is 
chosen with the above property then M exists. It then follows 
from equations (3.18) through (3.30) that A, must exist. We there- 
fore conclude from equations (3.31W) and (3.31V) that if A, = 0 
(or in practice less than some small constant multiplied by the cor- 
responding element of M) then Rq =0, but otherwise Ra =1. 

(Note that, as can be seen from equations (3.22) and (3.23), R^ is 
not used in any computations till A, has been determined) . Since 
the r-th diagonal element of R equals R^ (by definition of R'o ) , 
it follows that the determination of all matrices R^ in effect de- 
termines R. 


3 . 5 The Computation of the Incomplete Inverse in the Presence 

of Patterned Zeroes . 

Let us consider the symmetric matrix 



and its incomplete inverse 
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(3.47) 



D 


M* 




x; xi 


(3.48) 


Examination of elimination equations (3.18) through 
(3.26) show that A 2 and Y 2 do not modify Aj and . It can also 
be seen from the back-substitution equations (3.27) through (3.30) 
that Bj and X^ do not contribute anything to the computation of 
and X 2 • It follows that the partitions of M* (except for ) 
may equivalently be computed by the following . scheme : 


(i) 


Consider M, 


C 

Y. 


(3.49) 


and apply the elimination equations through A, obtaining 

C 




_T 

Y„ 


(3.50) 


(ii) 


Save Y 2 , A 2 and 


(iii) 


form M 

3 



(3.51) 


(iv) 


Obtain incomplete inverse of M, 


X, 


m: 




B 


1 J 


(3.52) 


(v) 


Save X j , B j and form 
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D 


(3.53) 


M, 


_T 


Y, A, 


(vi) Apply elimination equations through A, to obtain 


M 


5 




(3.54) 


The scheme as outlined above may also be applied when M 
is of the form 


C 

Y. 


M = 


0 

0 


0 

A, 

0 


0 

0 

'A, 


Y 

0 

0 

0 


(3.55) 


0 0 0 


Furthermore, the scheme may, in a fairly obvious way, be 
extended for the case of each A, (for i = 1, 2, . . . , n) itself being 

of the form (3.55). 


Note that the scheme just outlined is from a programming 

point of view far simpler (particularly when the matrix is stored in 

upper triangular form) than the conventional way of computing the 

true inverse (which is a special case of the incomplete inverse) 

through the scheme: 
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(i) 

Compute 

aV 

t 

a;^ 

Y 2 and 

C 

= C-Y 2 a;^ Y 2 . 

(3.56) 

(ii) 

Compute 

aV 

9 

aV 

Yj and 

D 

= (C-yJa;' Yj)“' 

(3.57) 

(iii) 

Compute 

Xi = 

-A’/ 

Y,D 

and Bj 

- 

A"' -X, (A"/ Y, y 

(3.58) 

(iv) 

Compute 

= 

-aV 

Y,D 

and B 2 

= 

-X 2 (a;' Y^y 

(3.59) 


When each A; itself is of the (3.55), the attractiveness of using 
the scheme based on equations (3.49) through (3.54) becomes even more 
apparent . 
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4.0 THE APPLICATION OF THE 


GENERALIZED INVERSE TO THE SOLUTION 
OF THE LINEAR LEAST SQUARES PROBLEM 


Let us consider the linear matrix equation 


u - AoX +e, 


(4.1) 


where 

the n-vector u represents a set of n observations, 

the m-vector x represents a set of m unknown parameters, 

the n-vector e represents the measurement noise of each 

observation, and 

the n m matrix Ao represents the mathematical modelling 

of the observations. 


The noise vector e is unknown. We do, however, assimie 
that its expected value (denoted by the operator E) vanishes, i.e., 

E(e) = 0 (4.2) 

It is also assumed that the covariance of the noise is 
known. Denoting this covariance by W”\ we thus have 

E(eeM - W~\ (4.3) 
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4.1 


The Solution 


It can be shown that if the measurement noise is random 
Gaussian then the most likely solution of equation (4.1) is obtained 
by minimizing the expression 

C(x) = (AoX-u/ W(AoX-u) (4.4) 

If X is the value of x, which minimizes C(x), then it 
follows that the partial derivations of C(x) with respect to x van- 
ish when X = X and 


AoW(AoX -u) = 0 


(4.5) 


Defining A and y, respectively, by 


A = aLwAo , 

and 


y = a'„Wu, 


it follows that 


A 

Ax = y, 


(4.6) 


(4.7) 


(4.8) 


and hence that 

X = A ^ y (4.9) 

X is known as the weighted least squares estimate. If the measure- 
ment noise is random Gaussian then it is also the maximum likelihood 


estimate. 
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4.2 The Error in the Solution . 

It follows from equation (4.1), (4.6), and (4. 

y = Ax + AgWe, 

where x represents the true parameter vector. From the 
equation (4.9) we obtain 


X -X = A~ ^ (aJ W e) 

We hence deduce from equation (4.2) that 

E(x) = X 


Since E( ee^ ) = W 


it follows from equation (^ 


E [(AlWe) (AoWe)'] = A. 

We hence deduce from equation (4.12) that 

E |^(x -x) (x -x/j = A"' , 

-I 

cov(x) = A 


4 . 3 The Value of C(x) . 

It follows from equations (4.4) and (4.5) that 
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.A, 

C(x) 


u^ Wu 


-u^ WA„ X 


) that 

(4.10) 
above and 

(4.11) 

(4.12) 
.6) that 


(4.12) 

(4.14) 



Defining Co by 


Co = u^Wu, 

it follows from equations (4.7), (4.9), and (4.14) that 


C(x) = Co -/A“V 


(4.15) 


(4.16) 


A . 4 Th e Computation of x, cov(x) , and C(x) Using the Incom - 

plete Inverse . 

Let the symmetric matrix M be defined by 


M = 


C„ 

y 


A 


(4.17) 


and R (corresponding to the same partitioning) by 



(4.18) 


It then follows from equations (2.59), (2.60), and (2.65) 


that 



(/A"' y -Co) 

A"^ y 



(4.19) 


From the above and equations (4.9), (4.13), and (4.16) we 
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Illlllll II I III 


II II 


I Ml 


■lll■■llll I 


II 


deduce that 

-C(x) x' 

A A (4.2) 

X cov(x) 

Note that -C(x), which is the upper left diagonal element 
of is obtained at the conclusion of the application of the elimi- 
nation algorithms. 



4. 5 The Solution When Some of the Parameter Values Already Are 

Known . 

Suppose Ao and x in equation (4.1) are partitioned 
according to 


Aq Aj , a.2 


(4.21) 


and 


X = 


(4.22) 


X 


so that equation (4.1) may be written in the form 


u = Aj Xj + A 2 X 2 + e 


(4.23) 


Suppose further that we already have an estimate X 2 for 
X 2 and that we do not wish to improve on that estimate. To reflect 
this, let us write equation (4.23) as 
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u' = AjXj + e' , 


(4.24) 


where 


u' = u -A^ , (4.25) 

.'incl 

e' = e - (x^ "^ 2 ^’ (4.26) 

The components of x^ are generally referred to as consid- 
er parameters , 

Corresponding to equations (4.6) through (4.9) the weighted 
least squares solution of equation (4.24) is given by 

x^ = (A\WAj)"' (AjWu’) (4.27) 

Corresponding to equation (4.11) we obtain for the error 
in the estimate 

Xj -X| = (a'^WAj) ' (A'^jWf') (4.28) 

From the above and equation (4.26) we obtain 

Xj -Xj = (a\wA^)"' a\ W e - a\ WA^ (x^ -x^ ) (4.29) 

We do not know the correct value of x^ . If we did, then 

A 

obviously we would have used that value for our estimate x^ . We 
do, however, assume that we know x^ to some known degree of accur- 
acy. Specifically, we assume that 

E(x^ ) = x^ , (4.30) 
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where is a known diagonal matrix. 

It follows from equations (4.29), (4.2), and (4.30) that 

E(x^) = x^ . (4.32) 

The sensitivity matrix Z is defined by 

Z = -(AjWAj)"' (a\wA^) (4.33) 

Equation (4.29) may then be expressed as 

Xj -X = (A^jWAj)"^ (A’^jWs) + Z(x^ (4.34) 

Since clearly E(ex 2 ) = 0 we deduce from equations (4.34), 
(4.3), and (4.31) that cov(x^), which is defined by 

cov(Xj) = E [ (xj -Xj)(Xi -Xj)’' j , (4.36) 

is given by 

cov(x^) = (aJwA^)"^ + (4.37) 

The Alias matrix L is defined by 

L = ZV^ 

Hence equation (4.36) may also be written as 
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cov(Xj) = (A'^WA^)"' + LL" (4.37) 

Let us now consider the more general form of equation 
(4.22) that generally occurs in practice. 

As before let us assume that we already have an estimate 
for some of the components of x, and also that we know the stand- 
ard deviation of the error in each estimate. 

Let S be defined as the diagonal matrix, whose diagonal 
elements equal either one or zero, such that an element equals one 
if the corresponding parameter is to be estimated but otherwise 
equals zero. 

If futher R is defined by 

R = I -S, (4.38) 

then we obtain corresponding to equation (4.23) 

u = AqRx + AoSx + e, (4.39) 

and corresponding to equations (4.24) through (4.26) 

u = AoR(x -x^) + e (4.40) 

where 

u = u - AoX,^ (4.41) 

and 

e = e + AoS(x -x^) , (4.42) 
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where is the initial or a priori estimate of x. 


Minimizing 


C(x) = [ ■ A„R(x-x^) J W - A„R(x-x^)j , (4.43) 

with respect to x, we obtain, when as before denoting the corres- 

A 

ponding value of x by x. 


Cw[ 


RAqW u - A^R(x-x. ) 


= 0 


(4.44) 


It follows from the above and equation (4.6) that 


j^y - AR(x-x^)j = 0‘ 


(4.45) 


where 


y = A^Wu 


(4.46) 


Since premultiplying equation (2.19) by R and then post- 
multiplying the resulting equation by R, yields, 


and 


RA*RA = R + RA*S, 


RA’^RAR = R 


(4.47) 

(4.48) 


we deduce after premultiplying equation (4.45) by RA* that 


R(x-x^) = RA*Ry 


(4.49) 


Since we already know Sx = (Sx^), equation (4.49) yields 
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the required least squares solution of equation (4.39). 


Let us next investigate the error in the estimate Rx. 
It follows from equations (4.46), (4.40), and (4.42) that 


y = a\w j^A„(x-x,,) + ej 


whence by equation (4.6), 


y = A(x-x,^) + A^We 


(4.50) 


From the above and equations (4.47) through (4.49) we de- 


duce that 


A 


R(x-x,^) = (R + RA*S)(x-x^) + RA*R A^We , 


and hence that 


R(x-x) = RA*R A^We - RA*S (x,, -x) (4.51) 

As before let us assume that 


E(Sx,,) = Sx , (4.52) 

and 

e[ S(x,,-x) (x,,-x)' S ] = SVVS , (4.53) 

where SV is a known diagonal matrix. 

It then follows that 
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Rx 


(4.54) 


A 

E (Rx) = 


and 

cov(Rx) = (RA*R)A(RA*R) + (RA* S)SWS (S/l’' R) 


The last equation may with the aid of equation (4.48) be 

written as 


cov(Rx) = R^ R + LL^ 
where the Alias matrix L is given by 


(4.55) 


L = -Rif SV (4.56) 

Finally we deduce from equations (4.43), (4.44), and 
(4.46) that 


C(x) = u^Wu - y^R(x-x^) 


(4.57) 


4.5.1 . The Computation of the Solution, The Alias and Covariance 

Matrices Using the Incomplete Inverse . 

Let the S3nnmetric matrix M be defined by 


M 


Co 

y 



(4.58) 


= u^ Wu . 


where- 
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C, 


(4.59) 



Let 1< and 5, when parti tioiied similarly to M 
tion (4.58) be defined by 



0 

0 


1 

o' 

R = 



, and S = 




0 

R 


0 

S 


Comparing the above with equations (3.1), (3.2) 
(3.10), (3.12), and (3.13) we deduce that 


where 


and 



X = (A*R-S)y 
D = -Co + y^ Rx 

B = A* 


Premultiplying equation (4.62) by R we deduce, 
aid of equation (449) that 


R(x-x^) = Rx , 

and hence from equations (4.63) and (4.57) that 

D = -C(x) 


in equa- 

(4.60) 
. (3.9), 

(4.61) 

(4.62) 

(4.63) 

(4.64) 
with the 

(4.65) 

(4.66) 


It follows from the above that 
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(4.67) 


The Alias matrix L and the covariance matrix cov(Rx) may 
be computed from equations (4.56) and (4.55), respectively. Note 
that since R and S are diagonal matrices, whose elements equal 
either zero or one, the elements of RA R and RA S are also elements 
of A*. Also note that it can easily be shown that the full covar- 
iance of X is given by 

cov(x) = RA*R + (SV+L)(SV+L)^ (4.68) 


4.5.2 The Change in the Solution When the Number of Parameters 
Being Estimated is Increased . 

In this subsection we shall investigate how the least 
squares solution of matrix equation (4.39) is changed, when one of 
the consider parameters is included in the set of parameters being 
estimated. 

Specifically, we shall consider how the solution is changed 
when S, as defined in the previous subsection, is changed to 
such that S-S^ is a diagonal matrix with only one non-zero element, 
which equals one. 

Corresponding to equation (4.38) let us define 

R^ = I-S' (4.69) 


and corresponding to equation(4 . 60) 

6o 


(4.70) 




0 

0 


0 


and 


1 0 
0 


It then follows from the definition of that 

r' = R+K , (4.71) 

where K is a diagonal matrix with only one non-zero element, which 
equals one. Let us assume that this is the (k+l)-th diagonal element, 
so that it is the k-th parameter that is now included in the set of 
parameters being estimated. 

Since u and, hence Co and y are independent of S, it 
follows that M as given by equation (4.58) is independent of S. 

However, M as given by equation (4.67) is the incomplete inverse of 

^ , 

M with respect to R, so that M will change as R changes from R . 

Retaining the notation M = M,R , we deduce, by comparison with equa- 
tions (2.144) through (2.146) and equation (2.155) that 

^ M*K + (I-K)j'^ [m*(I-K) -k] (4.72) 

Before futher considering equation (4.72), let us write M* 
in partitioned form as 


6i 



(4.73) 


M 


M 

m 

M. 


11 

T 

1 

'21 


m, 


U 




M, 

m, 

M 


22 


where the vector (nij , y , ) occupies the (k+l)th row of 

then follows that 


and 


M*K + (I-K) 


I nij 0 

0 u 0 

0 m2 ^ 


M*(I-K) -K = 


Mu 0 M2; 

m; -1 m2^ 

M 0 M 

1I21 ^‘■zz 


Inverting equation (4.74) we obtain 


^ M*K + (I-K) 




-y 


-m 


2 


0 

0 

I 


where 


and 
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Ej = mj/u 

E2 = mg/y 

y = -1/y 


It 


(4.74) 


(4.75) 


(4.76) 


(4.77) 

(4.78) 

(4.79) 



It follows from equations (4.72), (4.75), and (4.76) that 



-m^m^) 


(M2\ -mjm^) 



— T 

m, 

y 

— T 

m2 

(4.80) 

(Mji 

-m2mi^) 

^2 

(M 22 -mjm^) 



Equation (4.80) may be interpreted as follows. Let the 
estimate corresponding to be denoted by x' and let A*^i be de- 

noted by A*^. It then follows from equations (4.64), (4.67), (4.73), 
and (4.77) through (4.80) that 

(A*’),, = -1/(A*),, (4.81) 

(x'-x,,)^ = (x)k/(A*)kk (4.82) 

and A A A _ 

C(x’) = C(x) + (x'-x^)j, (x)^ (4.83) 

Of the three equations (4.81) through (4.83) the last one 
is of the greatest practical interest, because it tells us which para- 
meter, if included in the set of parameters to be estimated, yields 
the least value for C(x'). In other words, it tells us which addi- 
tional parameter to estimate in order to achieve the maximum improve- 
ment in the data fit. 


4. 6 The Solution When the Matrix is Computationally Singular . 

Two cases will be considered here: 

(i) Initial (or a priori) estimates are available for all 
parameters , and 


(ii) Initial (or a priori) estimates may be available for 
some parameters, but certainly not for all. 

In the next subsection, it will be shown how, in case (i) , 
the computational singularity may be avoided. 

4.6.1 Avoidance of the Computational Singularity When Independent 
A Priori Estimates are Available for All Parameters . 

Let us consider matrix equations (4.40) through (4.42). 
Since the a priori estimates x^ in a mathematical sense are equiva- 
lent to direct measurements of x, they are included in the vector u 
representing the observations. To reflect this, equations (4.39) 
through (4.42) may be written in partitioned form as 




Ao' 

Rx 

’aj" 

Sx 

'e' 


s 


+ 


+ 




R 


R 


1 

> 

1 

X 

1 


u' 


a; 

R(x-x,^ ) 

e' 


= 


+ 


0 


R 


R(x^ -x) 

L- ^ 


where 

u' = u' - AJx^ 

and 

e ' = e ' + Ao S (x-Xy^ ) 


(4.84) 


(4.85) 


(4.86) 


(4.87) 
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Comparing equations (4.85) and (4.40), we thus see that 



(4.88) 


u ' 


a: 


and Ao 


0 


R 


It also follows that the weighting matrix W is of the form. 


W 


W 0 
0 V"^ 


(4.89) 


where consistent with equations (4.52) and (4.53) 


E (x^ ) = X 


and 


e[(x^-x)(x^-x)^] = W 


(4.90) 


(4.91) 


From the above and equations (4.6) and (4.7) we deduce that 


A = A' + V~^ V"^ 


(4.92) 


where 


A’ = (Aj)' W'Ao’ 


(4.93) 


We similarly deduce from equation (4.46) that 


y = (Aj/ W'u' 


(4.94) 


Since by assumption the estimates of the components of x 


are independent it follows that V 


- 1 


is a diagonal matrix. 
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It was previously established (Subsection 3.4) that the 
matrix being inverted is computationally singular when = 0 in 
equation (3.19). Since the matrix must be at least semi-positive 
definite without any a priori estimates being available, A, should 
therefore be tested for being less than or equal to zero (or, in 
practice, some small predetermined constant. See (Morduch, 1975)) 
before the addition of the corresponding element of V”W~\ and if 
found to satisfy the test then both A^ and y, should be set to 
zero before the addition of the corresponding element of . 

yr should be set to zero, since if A, = 0 then the matrix cannot 
be semi-positive-definite unless y, = 0. 

Note that the procedure indicated above is mathematically 
equivalent to switching a parameter from being a 'solve for' to 
being a 'consider' parameter in order to be able to obtain a solu- 
tion. 


4.6.2 The Minimiim Norm Solution 

If it is found impossible to compute a definite solution 
to the least squares problem, then it is often desirable to obtain a 
minimiom norm solution. The computational procedure for a minimum 
norm solution is given in Subsection 2.6.1. In this subsection we 
shall discuss both the form of the norm and some justification for 
choosing a minimim norm solution. 

The norm of the estimate x, in accordance with equation 
(2.95) is given by 

A A A 

N^(x) = x^q'x 
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(4.95) 



I 


In as far as the components of x represent physical quan- 
tities, the terms of the sum appearing on the right hand side of 
equation (4.95) must have physically the same dimensions. It should 
mention that this condition is not, in general, met by the frequent- 
ly used Euclidean norm (Q=I) . A suitable norm is given by Q being 
the diagonal matrix, whose diagonal elements satisfy 

Qm = A„ (4.96) 

That this leads to a dimensionally correct norm is obvious since it 
has the same dimensions as C(x), the quantity being minimized for 
the least squares solution. It can be shown that the physical in- 
terpretation of the norm given by equation (4.96) is that not only 
is the weighted stun of the squares of the residuals minimized, but 
also the weighted sum of the squares of all individual contributions 
to the observation vector u. The proposed norm is quite general. 
Other norms that better fit particualr physical situations may, 
however, also be considered. 

Let us now discuss the justification for choosing a mini- 
mum norm solution. First, let us say that if the only requirement 
is that the weighted stun of the squares of the residuals should be 
a minimtun and no other information is given, then any solution of 
the least squares problem is as likely to be correct as any other. 
Such, however, is never the case when we are dealing with physical 
quantitities , which are always bounded. Also, by minimizing the 
corrections to the unknown parameter values, we ensure that, to the 
extent that they are in error, we also minimize any adverse effects 
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those parameters might have, when used to predict another set of 
observations . 


4. 7 The Expected Value of the Stim of the Squares of the 

Weighted Residuals . 

In this subsection we shall consider the expected value of 
C(x) as given by equation (4,57), viz . 


C(x) = u Wu -y R(x-x^) 


(4.57) 


It follows from the above and equations (4.49) and 
(4.46) that 


C(x) » u^W [u - AoRA*RA„Wu] 


(4.97) 


Since by equations (4.40) and (4.42) 


u = A, (x-x. ) + e , 


(4.98) 


we deduce with the aid of equations (4.6) and (4.47) that 

C(x) = u w|^Ao(x-x^) +e- Ao (RA* S+R) (x-x^ ) 

f^We] 


- A„RA"RA' 


i.e. , 
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C(x) = u'w[^A, (S-RA*S) (x-x^) + e -A, RA* RA"„ We] (4.99) 



and Q both PQ and QP 


Since, if for any two matrices P 
are defined, then trace(PQ) = trace (QP) , it follows that 

C(x) = trace |(x-x^ )u^ WAo (S-RA* S)| + trace jeti w| 

trace | a'„ Wcu^ WA„ RA r| (4.100) 

From the above and equations (4.98), (4.2), (4.3), (4.6), 
and (4.53) we obtain; 


E^C(x)^ = trace {wA(S-RA*S)j + trace |w“^ wj 


7 


trace |A RA*R 


(4.101) 


Post-multiplying equation (2.17) by S we deduce that 


AS - ARA*S = -SA*S 


(4.102) 


Since the diagonal elements of SA R vanish we find with 
the aid of equation (2.23) that 


trace 


j ara*r| 


= trace 


w 


(4.103) 


Hence , 


C^C(x)^ = trace | 


v-i 


W W 


trace 


|r| - trace jw(Si?S)| (4.104) 


If the number of observations equals n and the number of 
parameters being estimated is m, then 
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trace 


|w‘* wj 


= n and trace 


= m. (4.105) 


We hence obtain the result 


e ( c ( x )^ = n-m - trace jw(SA*S) 


(4.106) 
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5 . 0 SUMMARY 


The generalized matrix product of two square matrices A 
and B, denoted by A*B, is defined by 


A-B = -ARB + AS + SB 


( 2 . 16 ) 


where 


and 


R = R 


RR = R 


S = I-R 


( 2 . 2 ) 

( 2 . 3 ) 

( 2 . 12 ) 


It is shown that for any three square matrices A, B, and 


(A-B)-C = A-(B-C) 


( 2 . 4 ) 


and for any square matrix A 


A-(-R) = (-R)-A = A 


( 2 . 5 ) 


The generalized inverse of A with respect to R is denoted 
by A*R or, when no risk of confusion arises, simply by A*. It 
satisfies 


A -A* = A* - A = -R 


( 2 . 6 ) 
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Also 


(A*)* = A (2.8) 

and 

(A*)' = (aM* (2.11) 

When R defining the generalized inverse is diagonal, the 
inverse is referred to as the Incomplete inverse with respect to R. 

Given the matrix equation 

Ax = y (2.69) 

where A is a symmetric matrix and x and y are vectors, it is 
shown that there exists an R such that A (=A,R) is defined and 

SA*S = 0 (2.82) 

The general solution of equation (2.69) is given by 


X = Xo + (S-RA*'''S)a , (2.87) 

where 

Xo = RA*Ry , (2.86)' 

and the vector a satisfies Ra = 0, but is otherwise arbitrary. 
(Note that if A is non-singular then R = I, S = 0 and ^ ). 

The minimum norm solution Xq is obtained by minimizing the 

12 


norm 


Nq(x) = x’Q"‘x , (2.95) 

with respect to a. It is given by 

Xq = x„ + (S-H) (h"q-'H + Q-' )“'h"Q-'x„ , (2.115) 

or alternatively by 

Xq = (Q + QH')(HQH' +Q)"'x„ , (2.116) 


where 


H = RA*S 


A procedure for obtaining the incomplete inverse of a sym- 
metric matrix M is given in Section 3. It is shown that the pro- 
cedure may always be applied when (RMR+S) is positive-definite. 


Consider the linear matrix equation 


u = AoR(x-x^) + e 


where 


u = u -Ao 
e = e + Ao S (x-x^ ) 


(4.40) 


(4.41) 

(4.42) 
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u is a vector of observed measurements, e is the measurement 
noise, x is a parameter vector, is the a priori estimate of 

X, R and S are diagonal matrices satisfying equations (2,3) and 
(2.12), respectively. The components of x are either 'solve' or 
'consider' parameters. R and S are further defined such that 

Rii = 1 and S|, =0 if Xj is a 'solve' parameter, 

and 

Rii = 0 and Si, =1 if Xi is a 'consider' parameter. 


The weighted sum of the square of the residuals is given by 

C(x) = l^u -AoR(x-Xa)J W^u -AoR(x-x^)j (4.43) 


where 


W“' = E(ee" ) 


(4.3) 


C(x) is a minimum when x = x. x is the weighted least 
squares solution of equation (4.40). 


The symmetric matrix M is defined by 


M = 


Cc 

y 


— T 

y 


(4.58) 


where 


_T _ 

o = u Wu 
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C, 


(4.59) 



y = AqWu (4.46) 

and 

A = A^WA., (4.6) 


Corresponding to M, R and S are defined by 



1 

o 

o 
1 


o 

iH 
1 

R = 


and S = 



0 R 


[o sj 


(4.60) 


It is shown that 


M 


* 



(4,61), 

(4.64), 

(4.66) 


R(x-x^) = Rx , (4.65) 

and 

cov(Rx) = RA*R + LL^ , (4.55) 

where 

L = -RA*SV , (4.56) 

and 

SWS = e[s(x^-x) (x;,-x/ sj (4.53) 
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It is further shown that 

e[c(x)] = n-m - tracej VV(SA*S) I (4.106) 

where n is the number of observations and m is the number of 
'solve' parameters. Note that a priori estimates for 'solve' para- 
meters should be counted as observations 

If, e.g., is switched from being a 'consider' to a 

'solve' parameter, and the new estimate for x is denoted by x' , 
then 

C(x') = C(x) + (x); / (A*),, (4.83) 

(4.82) 
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